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Using a time quantified Monte Carlo scheme we performed simulations of the switching 
time distribution of single mono-domain particles in the Stoner-Wohlfarth approximation. 
We considered uniaxial anisotropy and different conditions for the external applied field. The 
results obtained show the switching time distribution can be well described by two relaxation 
times, either when the applied field is parallel to the easy axis or for an oblique external field 
and a larger damping constant. We found that in the low barrier limit these relaxation times 
are in very good agreement with analytical results obtained from solutions of the Fokker- 
Planck equation related to this problem. When the damping is small and the applied field 
is oblique the shape of the distribution curves shows several peaks and resonance effects. 

PACS numbers: 75.40. Mg, 75.50.Tt 



I. INTRODUCTION 

The study of magnetic nanoparticles is very interesting from the theoretical and experimental 
point of view, in addition to its important technological applications in magnetic data storage. 
Magnetic media for data storage are composed of tiny magnetic grains, like magnetic particles, 
which have to be nipped in the magnetic recording process. Increasing the storage density implies 
a reduction in the size of the magnetic grains. If the grains are too small they lose thermal stability 
reaching the so called superparamagnetic limit. The stability problem can be overcome by increasing 
the anisotropy of the grains but then higher fields (which are difficult to reach in current devices) 
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are needed to switch the magnetization in the recording process. In this sense a challenging issue 
about magnetization reversal arises which is how to have short reversal times, keeping the switching 
fields in small values [l[ 2] and controlling the effects of thermal fluctuations. The thermal stability 
problem of ferromagnetic particles was first studied by Neel and then by Brown 3, ^,[5] in the single 



spin approximation. Experimental studies for the switching time in a single magnetic nanoparticle 
Q, Q, Q] suggest in several cases the correctness of the Neel-Brown approximation. 

In his work Brown Q| derived the Fokker-Planck equation (FPE) for an assembly of parti- 
cles using the stochastic Landau-Lifshitz-Gilbert dynamics (sLLG) and then obtained analytical 
expressions for the lowest eigenvalue of the Sturm-Liouville problem associated with the FPE in 
asymptotic cases. The lowest eigenvalue of the FPE is associated to the longer relaxation time 
which in several regimes is supposed to control the switching process. Since then several works 
lave been done to obtain the lowest eigenvalue of FPE, with applied fields parallel js] and oblique 
id . lH|| to the easy axis. While in the past the attention was focused on the long time regime of 
the switching process, with the growing complexity of experimental devices and the need of ever 
smaller switching times, it turns to be desirable to have a knowledge of the whole process, from 
the earliest times to the longest ones, in order to explore alternative mechanisms for magnetization 
switching. 

Since the FPE cannot be solved analytically except in limiting cases, the short relaxation times in 



the early stages of the switching process has to be studied by numerical integration of the sLLG [12| 
equation, or solving numerically the FPE. For the intermediate time scales numerical integration 
of Langevin dynamics is not useful since the small time step needed in the numerical integration 
does not allow to extend the simulation beyond the range of nanoseconds. In order to overcome 



U, Ha, UM that 



this problem time quantified Monte Carlo methods (TQMC) were developed 
allow to extend the time span with a considerable gain in computational and programming efforts. 
In addition, Chubykalo et al. have also proved that these methods may be useful in the low 
damping condition, whenever the applied field is parallel to the anisotropy easy axis. However 
the method fails under an oblique applied field and low damping condition where precessional 



movement cannot be neglected. Recently, Cheng et al. [I7j mapped the sLLG dynamics and the 
Monte Carlo scheme through a Fokker-Planck approach introducing the precessional step in the 
TQMC simulations. With this improvement in principle TQMC simulations could be used also in 
the short times scales irrespectively of the configuration of the applied field and the value of the 
damping constant. However this improvement in the accuraccy is in detriment of the extent of time 
of the simulation. The situation is even worse when the system consists of a set of interaccting 
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particles IS}]. Nevertheless, even in these cases the TQMC scheme seems a better alternative than 
solving numerically the sLLG equation. The main advantage is a far simpler implementation, 
whithout the need to control the stability of the algorithm as when solving a differential equation 
numerically. The second advantage is that, even though the gain in the time span simulated relative 
to the sLLG equation is certainly not as large as in the strong damping regime where precession 
can be ignored, the TQMC scheme allows to adjust the time step in a less constrained manner than 
the sLLG equation, resulting in a real gain in the total times that can be reached by the simulation. 
The extent of this gain depends on the particular problem considered. 

In this work we used this approach to obtain the distribution of switching times in order to 
explore the incidence of short relaxation times in the switching process as a function of the damping 
and for different applied fields. We found that for external fields parallel to the easy axis the Brown 
approach works very well in practically the whole time span. The numerical results are extremely 
well described with only the two largest relaxation times from the solution of the Fokker-Planck 
equation. When the external field is oblique there are no analytical solutions available, except for 
the largest time scale. In this case, we found that two relaxation times are enough to describe the 
distribution of times when the damping is high. When the damping is low, the switching mechanism 
is dominated by precession, and the short time behavior is more complex, showing several peaks, 
which reflect the presence of different resonance frequencies. 



II. THEORETICAL BACKGROUND 



A single-domain particle can be modeled in the Stoner-Wohlfarth (SW) approximation where 
the magnetic state of the particle is described by a single magnetic moment fh whose strength is 
equal to the total magnetic moment of the particle \m\ = M s v. Here M s is the magnetization of 
saturation of the particle and v is the particle volume. In the SW model the energy density V of a 
particle with uniaxial anisotropy under an external applied field is expressed as: 



-PV = a[(n-s) 2 + 2h-s\, (1) 

here j3 = v/ksT, s = M/M s and n are unit vectors defining the magnetization and the easy 
axis direction, respectively The applied field h is expressed in reduced units h = H/Hk, with 
iffe = 2K/hqM s being the anisotropy field, a = Kv/ksT is a dimensionless constant where K is 
the anisotropy constant. In a classical approximation, the dynamics of a reduced magnetic moment 



4 



s under thermal fluctuations is modeled by the stochastic Landau-Lifshitz-Gilbert (sLLG) equation, 

§ = TT^ * x + ^ ~ a * x + ^ )] ' (2) 

where 7 is the gyromagnetic ratio, a = rp/M s is a dimensionless damping coefficient and rj is the 
damping coefficient in Gilbert's equation. The effective field h e ff is given by the particle energy 
gradient, 



- (3) 

The stochastic fluctuating field hfi is assumed as a Gaussian stochastic process with the following 
statistical properties: 



<h f i ;i >=0, < hfi ti (t)hfi tj (s) >= (j,6 itj 6(t- s), (4) 

where i and j stand for the cartesian components. 

In spherical coordinates, the FPE for the probability W(6, (f>, t) of finding one the magnetic 
moment at time t within a solid angle dO is given by [4, 3|: 



dW 1 „ 2rxr m2irm b fdVdW 8VdW\ 

-V 2 W + abV 2 V W + — — — — (5) 



dt 2t n sm{9) V d6 d<p d</> 89 

+a \d6~d0 + sin(#) 2 + ~d0~d0 

where the Neel time r^ 1 = n^ 2 (l + a 2 ) is a characteristic diffusional time and b = j/(l + a 2 )M s . On 
the other hand, in order to satisfy the equilibrium statistical properties in the stationary regime, 



(1 + a 2 ) 1 f(t , 
T N = a — -. 6 

In general solutions of © can not be found analytically. However, the relaxation of any initial 
probability state can be formally described by a sum of exponentials: 



00 

W(9, (t>) = W + Y, A n Fn(0, <j>) eM-t/n), 

n=l 



(7) 



where Tj are related to the eigenvalues of the Sturm-Liouville associated problem 
to 



according 



2t n 
Xi ' 



(8) 



Besides the diffusional Neel time, the other characteristic time scales in this problem are the pre- 
cessional time r v = {^^lqHk/{^ + a 2 )) -1 and the damping time tk = T p /a |19|. Then, eq. ([8]) 
expressed in terms of the damping times becomes t^tk] = 



III. MONTE CARLO SIMULATIONS 



In our simulations we use a hybrid Monte Carlo method 
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171 ] that emulates the stochastic 



dynamics of the LLG equation. This method combines the Metropolis or heat bath MC scheme 



with a random displacement of the magnetic moment within a cone [13] and a precessional spin 
motion. The random displacement is obtained by adding to the normalized magnetic moment a 
random vector uniformly distributed within an sphere of radius R <C 1 and then normalizing the 
resulting vector again, while the magnitude and direction of the precessional motion is given by 



As = -$sx h e ff, (9) 

where <I> = j^R 2 . In this MC scheme the magnetic moment update is chosen with equal probability 
between a precessional step and a random displacement. The acceptance rate A(f3AV) of the 
random motion is based in the heat bath procedure, 



A(0AV) = 1/[1 + exp(/3AV)]. (10) 

By means of a detailed comparison between the Fokker-Planck equation representing the MC 
stochastic dynamics and the corresponding equation associated with the LLG micromagnetic equa- 
tion, it is possible to obtain a very accurate mapping between Monte Carlo Steps and the real time 
scale from the LLG equation jl7] through the following relation: 



At[r K ] = —R 2 At[MCS], 
where the real time scale is expressed in units of the damping characteristic time tr. 



(11) 
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IV. RESULTS 

We have performed numerical evaluations of the switching time distribution with the following 
protocol: we start with the magnetic moment pointing in the easy axis direction, in our case 
the z axis, and then we apply an inverse magnetic field of different strengths and direction. The 
switching time is defined as the time required for the z component of the magnetic moment to 
change its sign. Since we want to compare our results with analytic predictions from the Fokker- 
Planck equation [5j we count every time the magnetic moment crosses the equatorial line, i.e., we 
compute the whole probability that the magnetic moment attains an angle 9 = tt/2 in a given time 
t. Otherwise we would be computing the first passage time, for which there are still less analytic 
results available. We performed 10 6 realizations to obtain the switching time distribution P(t). 
bmce P(t) oc Jq W{6 = tt/2, <j>, t) d(f>, it has the same relaxational behavior than W(6,t). 

In order to test the confidence of our results we simulated switching times distribution for the 
low barrier limit (a < 1), with the applied field parallel to the easy axis. In this case the two 
smallest eigenvalues are given by 



Ai =2-|(l -h 2 a)a + J|a 2 + tf(a 3 ) (12) 
5 o75 

A 2 = 6 - ~(1 - h 2 a)a + —a 2 + tf(a 3 ). 
7 343 

From equations (JSj) and ffTTj) . and considering that a time of 3t\ is enough to obtain a complete 
distribution curve, the number of MCS that should be used for each realization is: 



At[MCS] = -JjL. (13) 

In the simulations we used R = 0.03, which is a good compromise between accuracy and efficiency 
In figured] we present the switching time distribution for an external field h = 0.292 parallel to the 
easy axis and low energy barrier a = 1. At short times the probability that the particle switches 
is zero since there is a minimal time necessary to surmount the barrier. Except for the very short 
times, the distribution is well fitted using two exponentials P(t) = a^+ai exp(— tjr\)-\-ai exp(— t/T-z) 
where the relaxation times t\ and r 2 are obtained through eq. J8)) using the eigenvalues given in eqs. 
(fT2l) . We can see that two relaxation times are enough to obtain a good fit of the switching time 
distribution in a wide range and a very good agreement for the relaxation times obtained through 
the FPE. The switching time distribution at large times is finite because the energy barrier is small 
and the particle attains thermodynamic equilibrium with a finite probability of being at 6 = tt/2. 
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MC simulation 
a + a 1 exp(-t/t 1 ) - a 2 exp(-t/t 2 ) 
a - a 2 exp(-t/t 2 ) 
a., exp(-t/t 1 




FIG. 1: Switching time distribution in the low barrier case a — 1 and parallel applied field h = 0.292, here 
the damping constant is a = 1. The relaxation times used in the fitting r\ = 7.71 and T2 = 1.79 are obtained 
from eqs. (O) and lfl2|) . 



In fig. [21 the energy barrier is increased (a = 10) keeping the field parallel, with h = 0.4. In this 
case the switching time distribution has a peak before relaxing to zero. In this case the distribution 
goes to zero at large times since the magnetic moment gets trapped in the deepest minimum A 



similar behavior is found by integrating the sLLG equation 



12). Like in fig. [H this curve is well 



fitted by two relaxation times. The longer time t\ used in this fitting corresponds to results of 



Coffey et al. 



III]. The secondary relaxation time T2 , which is supposed to be related to the second 



eigenvalue is much lower than t\ and is important in the very first stages of the relaxation, having 
little influence in the switching mean time < t >. We also show in fig. [2] two curves corresponding 
to small and large damping constants a = 0.1, 100. These curves are indistinguishable within the 
statistical error, when plotted in units of the damping time tk- This is because if the applied field 
is parallel to the easy axis, then the potential energy has azimuthal symmetry and the precessional 
motion, which is important for small damping, has no influence in crossing the barrier. We will 
see below that a different situation is observed if the applied field is oblique. However, from the 
point of view of the Monte Carlo simulation, the damping constant has a critical influence since 
the precessional factor <3? has to be kept at small values in order to correctly follow the trajectory. 
Then, reducing the damping constant implies a reduction of R in eq. j9|), and the number of MCS 
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FIG. 2: High barrier switching time distribution with a = 10 for two damping constants a = 0.1 and 100 
and parallel applied field h = 0.4. The longer relaxation time n = 52.0 is obtained from table I [ll|. The 
secondary relaxation time is obtained from fitting: t-i = 1.72 and the mean time t m — 42.5. 

required in the simulation notably increase (see eq. (fl~3l) ). 

Figure [3] is similar to fig. (2J h = 0.4, but now the applied field is at an angle of 7r/4 with respect 
to easy axis. The damping constant has now influence on the switching behavior, this is due to the 
coupling between the longitudinal and normal modes [^J. When the applied field is oblique the 
particle crosses the barrier through a saddle point and the precessional motion influences the way 
the energy landscape is explored. If the damping constant is high the particle tends to relax to the 
nearest minimum and only through thermal fluctuations the saddle point can be found, wheres if 
the damping constant is kept small, precessional motion affects the switching process. Note that a 
shoulder is present after the first peak in the case of small damping. 

In fig. 0]the damping constant is decreased even more to a = 0.01. Now new peaks appear in the 
switching time distribution, which shows resonance-like effects. If the starting energy is similar to 
the energy of the saddle point and the damping constant is low enough, the particle does not relax 
quickly and keeps its energy nearly constant for a long time, letting the magnetic moment to cross 
the equatorial line more than one time. This behavior is analyzed in more detail in figure EJ where 
the curve of fig. [4] is plotted together with the distribution of the first and second passage times. 
The figure shows that the first and second peaks in the switching time distribution correspond to 
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FIG. 3: High energy barrier a = 10 under an oblique applied field <j) = ir/4, h — 0.4 for two different values 
of the damping constant a = 0.1, 100. It can be seen that in the larger damping case the distribution is well 
described by two relaxation times whereas in the low damping case a secondary peak appears. 



the main peaks of the first and second passage times distributions, respectively. These distributions 
show also secondary peaks, which probably correspond to different paths of switching in the energy 
landscape. From the figure is also clear that the first passage time probability goes to zero at 
t/TK ~ 70, while the whole distribution stays finite until much longer times. This fact means that 
the magnetic moment can go back and forth across the saddle point and the magnetic moment 
keeps switching between the basins of the two minima during a long time. Although the barrier is 
high, the small damping makes the magnetic moment to follow a long trajectory before settling in 
the final state. 



V. SUMMARY 



In summary, by performing Monte Carlo simulations including a precessional step, we have 
obtained switching time distributions of magnetic particles in the Stoner-Wolfarth limit, for different 
configurations of the applied field, different values of the damping constant and different heights 
of the energy barriers. We conclude that if the damping constant is high enough the distributions 
are well described by two relaxation times associated with the eigenvalues of the Fokker-Planck 
solutions of the corresponding Landau-Lifshitz-Gilbert dynamics, irrespectively of the configuration 
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FIG. 4: Switching time distribution for a very low damping value a = 0.01. Like in fig. [3] here a = 10, 
(j> — 7r/4 and h = 0.4. We can see at least four peaks and an exponential relaxing behavior at the longest 
times. 
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FIG. 5: Switching time distribution of fig. |4] together with the first and second passage time distribu- 
tion. Clearly the first and second peak of the total distribution correspond to the first a second passage, 
respectively. 
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of the applied field. If the damping constant takes small values and the applied field is oblique to 
the anisotropy axis, the distributions show resonance effects, evidencing the importance of the 
precessional motion in the inversion mechanism. The present Monte Carlo algorithm allows to 
study these precessional effects in detail, without the need to solve the LLG equations. We showed 
that the first two peaks in the distribution functions correspond to the first and second passage 
times of the magnetization through the equator. In all cases the characteristic inversion time is 
given by the smallest eigenvalue of the FPE, while in the cases with strong damping the whole 
distribution can be very well described by only two relaxation times. 
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